﻿using System;
using System.Linq;

namespace MathMethods
{
    /// <summary>
    /// Методы Решения СЛАУ
    /// </summary>
    public static class NumericalMethods
    {
        /// <summary>
        /// Метод Гаусса с выбором главного элемента (прямой метод).
        /// </summary>
        /// <param name="mm">Матрица, задающая СЛАУ класса Matrix c вектором правых частей.</param>
        /// <returns>Вектор решений.</returns>
        public static Row Gauss(Matrix mm)
        {
            /// Масска для восстановления изначального порядка
            /// неизвестных
            int[] mask = new int[mm.rowCount];
            for (int i = 0; i < mm.rowCount; i++)
                mask[i] = i;

            /// Прямой проход
            for (int i = 0; i < mm.rowCount; i++)
            {
                int maxIndx = i;
                double maxval = mm[i, i];

                /// Поиск главного элемента в строке
                for (int j = i + 1; j < mm.colCount - 1; j++)
                    if (Math.Abs(maxval) < Math.Abs(mm[i, j]))
                    {
                        maxval = mm[i, j];
                        maxIndx = j;
                    }

                if (maxval == 0)
                {
                    Row x = new Row(mm.rowCount);
                    for (int ii = 0; ii < mm.rowCount; ii++)
                        x[ii] = double.NaN;
                    return x;
                }

                /// Перестановка столбцов
                if (i != maxIndx)
                {
                    mm.ExchangeColumns(i, maxIndx);
                    int tmp = mask[i];
                    mask[i] = mask[maxIndx];
                    mask[maxIndx] = tmp;
                }

                /// Делим строчку на главный элемент
                mm[i] = mm[i] / maxval;

                /// Вычетаем текущую строчку
                /// (предварительно домножив) из ниже стоящих строк
                for (int j = i + 1; j < mm.rowCount; j++)
                    mm[j] = mm[j] - mm[j, i] * mm[i];
            }

            // Обратный проход
            for (int i = mm.rowCount - 1; i >= 0; i--)
            {
                /// Вычетаем текущую строчку
                /// (предварительно домножив) из выше стоящих строк
                for (int j = i - 1; j >= 0; j--)
                    mm[j] = mm[j] - mm[j, i] * mm[i];
            }

            /// Восстановление порядка вхождения неизвестных
            /// и формирование результирующего вектора
            Row res = new Row(mm.rowCount);
            for (int i = 0; i < mm.rowCount; i++)
                res[mask[i]] = mm[i, mm.colCount - 1];
            return res;
        }

        /// <summary>
        /// Поиск решения СЛАУ методом Зейделя (итерационный метод).
        /// </summary>
        /// <param name="mm">Матрица, задающая СЛАУ класса Matrix c вектором правых частей.</param>
        /// <param name="Eps">Точность вычисления.</param>
        /// <returns>Вектор решений.</returns>
        public static Row Seidel(Matrix mm, double Eps)
        {
            Row x = new Row(mm.rowCount);
            for (int i = 0; i < mm.rowCount; i++)
                x[i] = mm[i, mm.colCount - 1];

            Row prev_x = new Row(mm.rowCount);
            bool tmp;
            int counter = 0;
            do
            {
                counter++;
                if (counter > 200)
                {
                    for (int i = 0; i < mm.rowCount; i++)
                        x[i] = double.NaN;
                    Console.WriteLine("Произошло зацикливание");
                    return x;
                }
                tmp = true;
                for (int i = 0; i < mm.rowCount; i++)
                {
                    double var = 0;
                    for (int j = 0; j < mm.rowCount; j++)
                        if (i != j)
                            var += mm[i, j] * x[j];
                    prev_x[i] = x[i];

                    x[i] = (mm[i, mm.colCount - 1] - var) / mm[i, i];
                    tmp = tmp && (Math.Abs(prev_x[i] - x[i]) < Eps);
                }
            } while (!tmp);

            return x;
        }

        /// <summary>
        /// Вычисление определителя квадратной матрицы Matrix.
        /// </summary>
        /// <param name="mm">Матрица, задающая СЛАУ класса Matrix c вектором правых частей.</param>
        /// <returns>Определитель матрицы.</returns>
        public static double SquareMatrixDeterminant(Matrix mm)
        {
            double det = 1;
            for (int i = 0; i < mm.rowCount; i++)
            {
                det *= mm[i, i];
                mm[i] /= mm[i, i];
                for (int j = i + 1; j < mm.rowCount; j++)
                {
                    det *= mm[j, i];
                    mm[j] /= mm[j, i];
                    mm[j] -= mm[i];
                }
            }
            det *= mm[mm.rowCount - 1, mm.rowCount - 1];
            return det;
        }

        /// <summary>
        /// Метод Кремера для решения СЛАУ.
        /// </summary>
        /// <param name="mm">Матрица, задающая СЛАУ класса Matrix c вектором правых частей.</param>
        /// <returns>Вектор решений.</returns>
        public static Row CramersRule(Matrix mm)
        {
            Row x = new Row(mm.rowCount);
            Matrix temp = new Matrix(mm.rowCount, mm.rowCount);
            for (int i = 0; i < mm.rowCount; i++)
                for (int j = 0; j < mm.rowCount; j++)
                    temp[i, j] = mm[i, j];
            double det = SquareMatrixDeterminant(temp);

            for (int k = 0; k < mm.rowCount; k++) // который корень считаем
            {
                for (int i = 0; i < mm.rowCount; i++) // i-тая строка
                    for (int j = 0; j < mm.rowCount; j++) // j-тый столбец
                    {
                        if (k == j)
                            temp[i, j] = mm[i, mm.colCount - 1];
                        else
                            temp[i, j] = mm[i, j];
                    }
                x[k] = SquareMatrixDeterminant(temp) / det;
            }
            return x;
        }

        /// <summary>
        /// Прямой степенной метод.
        /// </summary>
        /// <param name="mm">Симметричная матрица.</param>
        /// <param name="initialVector">Вектор начальных приближений.</param>
        /// <param name="epsilon">Точность вычислений.</param>
        /// <returns>Спектральный радиус матрицы.</returns>
        public static double PowerMethod(Matrix mm, Matrix initialVector, double epsilon)
        {
            // начальное приближение
            Matrix y_old = initialVector;
            Matrix y_new = MultiplyMatrix(mm, y_old);
            int i = 0;
            double lambda1 = 0, lambda2 = 0;
            do
            {
                i++;
                lambda1 = NormVector(y_new) / NormVector(y_old);
                y_old = y_new;
                y_new = MultiplyMatrix(mm, y_old);
                lambda2 = NormVector(y_new) / NormVector(y_old);
            }
            while (Math.Abs(lambda1 - lambda2) > epsilon);
            return lambda2;
        }

        /// <summary>
        /// Интерполяция методом Лагранжа.
        /// </summary>
        /// <param name="x">В этой точке мы ищем значение функции.</param>
        /// <param name="xi">Массив известных точек.</param>
        /// <param name="yi">Массив известных значений функции в известных точках.</param>
        /// <returns>Значение функции в точке х.</returns>
        public static double LagrangePolynomial(double x, double[] xi, double[] yi)
        {
            if (xi.Length != yi.Length)
                return double.NaN;

            foreach (double temp in xi)
                if (Array.IndexOf(xi, temp) != Array.LastIndexOf(xi, temp))
                    return double.NaN;

            int length = xi.Length;
            double res = 0;

            for (int i = 0; i < length; i++)
            {
                double basisPolynomial = 1;
                for (int j = 0; j < length; j++)
                    if (i != j)
                        basisPolynomial *= (x - xi[j]) / (xi[i] - xi[j]);
                res += yi[i] * basisPolynomial;
            }
            return res;
        }

        /// <summary>
        /// Вычисление нормы вектора как максимальное по модулю число.
        /// </summary>
        /// <param name="mm">Метрица-вектор.</param>
        /// <returns>Норму.</returns>
        public static double NormVector(Matrix mm)
        {
            double norm = Math.Abs(mm[0, 0]);
            if (mm.colCount == 1)
            {
                for (int i = 0; i < mm.rowCount; i++)
                    if (Math.Abs(mm[i, 0]) > norm)
                        norm = Math.Abs(mm[i, 0]);
            }
            else
                for (int i = 0; i < mm.colCount; i++)
                    if (Math.Abs(mm[0, i]) > norm)
                        norm = Math.Abs(mm[0, i]);
            return norm;
        }

        /// <summary>
        /// Умножение двух матриц.
        /// </summary>
        /// <param name="mm">Мервая матрица MxN</param>
        /// <param name="ll">Вторая матрица NxQ</param>
        /// <returns>Матрица MxQ</returns>
        public static Matrix MultiplyMatrix(Matrix mm, Matrix ll)
        {
            Matrix res = new Matrix(ll.colCount, mm.rowCount);
            for (int i = 0; i < mm.rowCount; i++)
                for (int j = 0; j < ll.colCount; j++)
                {
                    double currentValue = 0;
                    for (int r = 0; r < mm.colCount; r++)
                        currentValue += mm[i, r] * ll[r, j];
                    res[i, j] = currentValue;
                }
            return res;
        }

        /// <summary>
        /// Вывод матрицы класса Matrix.
        /// </summary>
        /// <param name="mm">Матрица, задающая СЛАУ класса Matrix c вектором правых частей.</param>
        public static void PrintMatrix(Matrix mm)
        {
            Console.WriteLine("Matrix");
            for (int i = 0; i < mm.rowCount; i++)
            {
                for (int j = 0; j < mm.colCount; j++)
                    Console.Write(mm[i, j] + "  ");
                Console.WriteLine();
            }
        }

        /// <summary>
        /// Метод прогонки (он же алгоритм Томаса, он же tridiagonal matrix algorithm (TDMA))
        /// </summary>
        /// <param name="a">Нижняя диагональ</param>
        /// <param name="b">Главная диагональ</param>
        /// <param name="c">Верхняя диагональ</param>
        /// <param name="d">Вектор правых частей</param>
        public static void TDMA(double[] a, double[] b, double[] c, double[] d)
        {
            int size = a.Length;

            double[] p = new double[size + 1];
            double[] q = new double[size + 1];
            double[] x = new double[size + 1];

            p[1] = -c[0] / b[0];
            q[1] = d[0] / b[0];

            for (int i = 1; i < a.Length; i++)
            {
                p[i + 1] = -c[i] / (b[i] + a[i] * p[i]);
                q[i + 1] = (d[i] - a[i] * q[i]) / (b[i] + a[i] * p[i]);
            }

            x[size - 1] = q[size - 1];
            for (int i = 19; i > 0; i--)
                x[i - 1] = p[i] * x[i] + q[i];

            Console.WriteLine("Ответ:");
            for (int i = 0; i < a.Length; i++)
                Console.WriteLine("{0:f2}", x[i]);
        }

        /// <summary>
        /// Второй полином Ньютона для инторполяции в конце.
        /// </summary>
        /// <param name="x">Точка, значение в которой ищем.</param>
        /// <param name="xi">Массив аргументов.</param>
        /// <param name="yi">Массив значений функции.</param>
        /// <returns>Значение функции в точке.</returns>
        public static double NewtonPolynomial(double x, double[] xi, double[] yi)
        {
            // find step, it is always the same 
            // e.g. 0, 0.1, 0.2,..
            // here step is 0.1
            double h = xi[1] - xi[0];
            // just t
            double t = (x - xi[xi.Length - 1]) / h;
            int l = yi.Length - 2;
            //int g = 1;
            // start
            double result = yi[yi.Length - 1];
            for (int i = l, g = 1; i >= 0; i--, g++)
            {
                double koeff = 1;
                for (int w = 0; w < g; w++)
                    koeff *= t + (g - 1);
                koeff /= NumericalMethods.Factorial(g);
                result += koeff * NumericalMethods.FiniteDifference(g, yi);
            }
            return result;
        }

        /// <summary>
        /// Функция для полинома Ньютона.
        /// </summary>
        /// <param name="x">Аргумент.</param>
        /// <returns></returns>
        public static double FunctionForNewtonPolynomial(double x)
        {
            // переопределить для себя
            return Math.Log(1 + x * x);
        }

        /// <summary>
        /// Вычисление конечной разности.
        /// </summary>
        /// <param name="degree">Порядок k.</param>
        /// <param name="y">Значения функции.</param>
        /// <returns>Конечную разность.</returns>
        public static double FiniteDifference(int degree, double[] y)
        {
            double result = 0, koeff;
            int i = y.Length - 1 - degree, k = degree, g = 0;
            bool f = true;
            // start
            result += y[i + k];
            for (int q = i + k - 2; q >= i - 1; q--)
            {
                koeff = f ? -1 : 1;
                f = f ? false : true;
                if (q != i - 1)
                {
                    for (int w = 0; w < g + 1; w++)
                        if ((k - g) != 0)
                            koeff *= (k - g);
                    koeff /= NumericalMethods.Factorial(g + 1);
                }
                g++;
                result += koeff * y[q + 1];
            }
            return result;
        }

        /// <summary>
        /// Вычисление факториала числа.
        /// </summary>
        /// <param name="n">Число, факториал которого надо вычислить.</param>
        /// <returns>Значение факториала.</returns>
        public static long Factorial(int n)
        {
            long result = 1;
            for (int i = 2; i <= n; i++)
                result *= i;
            return result;
        }

        /// <summary>
        /// Выводит значение функции и ее аргумент в читаемом виде: y(x) = y.
        /// </summary>
        /// <param name="xi">Массив аргументов.</param>
        /// <param name="yi">Массив значений функции.</param>
        /// <param name="precision">Точность, сколько знаком после запятой.</param>
        public static void PrintFunctionValues(double[] xi, double[] yi, int precision)
        {
            if (xi.Length != yi.Length)
            {
                Console.WriteLine("Размерности массивов не совпадают!");
                return;
            }
            for (int i = 0; i < xi.Length; i++)
                Console.WriteLine("y(" + xi[i] + ") = " + yi[i].ToString("F" + precision + ""));
        }

        /// <summary>
        /// Интерполяционный метод наименьших квадратов.
        /// </summary>
        /// <param name="xi">Вектор аргументов.</param>
        /// <param name="yi">Вектор значений функции.</param>
        /// <param name="x">Аргумент.</param>
        /// <returns>Значение фугкции в искомой точке.</returns>
        public static double LeastSquares(double[] xi, double[] yi, double x)
        {
            if (xi.Length != yi.Length)
                throw new ArgumentException("Vector lengths do not match.");
            int n = xi.Length;
            double
                // вычисление формулы 1
                s1 = xi.Aggregate(0.0, (current, item) => current + item),
                // вычисление формулы 2
                s2 = yi.Aggregate(0.0, (current, item) => current + item),
                // вычисление формулы 3
                s3 = xi.Aggregate(0.0, (current, item) => current + item * item),
                s4 = 0;
            // вычисление формулы 4
            for (int i = 0; i < n; i++)
                s4 += xi[i] * yi[i];
            // вычисление коэффициентов прямой
            double
                a = (s4 - (s1 * s2) / n) / (s3 - s1 * s1 / n),
                b = (s2 - a * s1) / n;
            return a * x + b;
        }

        /// <summary>
        /// Stirling interpolation formula.
        /// </summary>
        /// <param name="x">Argument.</param>
        /// <param name="xi">Array of X arguments.</param>
        /// <param name="yi">Array of Y values, f(x) = y.</param>
        /// <returns>F(x).</returns>
        public static double StirlingInterpolation(double x, double[] xi, double[] yi)
        {
            // result
            double result = 0.0;
            // if lengths of x and y are not equal or
            // x and y lenght are not odd
            // return NaN result
            if (xi.Length != yi.Length || xi.Length % 2 == 0 || yi.Length % 2 == 0)
                return double.NaN;
            // check if step if constant
            if (!ConstantStep(xi))
                return double.NaN;
            // calculate step
            double step = xi[1] - xi[0];
            // center of array
            int center = (int)xi.Length / 2;
            // calculate q
            double q = (x - xi[center]) / step;
            for (int i = 1; i < xi.Length; i++)
            {
                double
                    // current function value
                    member = 0.0,
                    // temp value of q
                    temp_q = 1.0;
                // end value for q calculations
                int end = (int)Math.Ceiling((double)i / 2);
                // caculate q
                if (i > 2)
                {
                    for (int j = 1; j < end; j++)
                        temp_q *= q * q - j * j;
                }
                // calculate different members for different count
                if (i % 2 == 1)
                {
                    double temp_y = FiniteDifference(i, -end, yi) - FiniteDifference(i, -end + 1, yi);
                    member = (q * temp_q) * temp_y / (NumericalMethods.Factorial(i) * 2);
                }
                else
                {
                    double temp_y = FiniteDifference(i, -end, yi);
                    member = (q * q * temp_q) * temp_y / NumericalMethods.Factorial(i);
                }
                // add member to the result
                result += member;
            }
            // add y0 to the result
            result += yi[center];
            return result;
        }

        /// <summary>
        /// Counts finite difference.
        /// </summary>
        /// <param name="degree">Finite difference degree.</param>
        /// <param name="startIndex">Y index (usually is 0). But works for negative values too.</param>
        /// <param name="y">Initial array.</param>
        /// <returns>Finite difference.</returns>
        public static double FiniteDifference(int degree, int startIndex, double[] y)
        {
            // result
            double result = 0.0;
            // if start index is less then 0
            // and y.Length is odd,
            // we must find center index
            int center = 0;
            if (startIndex < 0 && y.Length % 2 == 1)
                center = (int)y.Length / 2;
            // get coefficients of finite difference
            int[] c = PascalTriangle(degree);
            // multiply odd coeffs by (-1)
            for (int i = 1; i < c.Length; i += 2)
                c[i] = -c[i];
            // calculate polinomial
            for (int i = degree + startIndex, coeff = 0; i >= startIndex; i--, coeff++)
                result += c[coeff] * y[i + center];
            // return result
            return result;
        }

        /// <summary>
        /// Generates coefficients for Pascal triangle.
        /// </summary>
        /// <param name="degree">Row number.</param>
        /// <returns>Coefficients.</returns>
        public static int[] PascalTriangle(int degree)
        {
            // coeffs array
            int[] c = new int[degree + 1];
            int i;
            // init all with 0
            for (i = 1; i <= degree; i++)
                c[i] = 0;
            c[0] = 1;
            // generate coeffs
            for (int j = 1; j <= degree; j++)
                for (i = j; i >= 1; i--)
                    c[i] = c[i - 1] + c[i];
            return c;
        }

        /// <summary>
        /// Chech input array if it has constant step.
        /// </summary>
        /// <param name="x">Array to check.</param>
        /// <returns>Bool true or false.</returns>
        public static bool ConstantStep(double[] x)
        {
            // init with true
            bool r = true;
            // get step
            double step = x[1] - x[0];
            // loop
            for (int i = 1; i < x.Length; i++)
            {
                if (x[i] - x[i - 1] != step)
                {
                    r = false;
                    break;
                }
            }
            return r;
        }
    }
}

